global name "s3_simulation"

global datafile "datacvr" 


global main "C:/Users/silvio/Documents/CVR/ryp"
global output "$main/output"  
global data "$main/data"  
global code "$main/code"  
global excel "$main/excel"  
global odata "$data/Intermuestra1-2"  

global coord "$data/coordinates"  
global geocoord "$data/geocoord"  


*1 Fix DGP. Simulation with parameters and best model from best data for SL
cd $main
capture log close
log using "$output/${name}.log", replace

matrix drop _all
set more off

set matsize 1392
scalar nstrata=58                      /*=58, CVR strata*/

set seed 98034

use "${data}/${datafile}" , replace

mkmat y
scalar perpeh=1

**************************************************************

scalar nsimu=1000
matrix one7=J(1,7,1)
matrix res=J(9, 16, 0)
******************************

scalar icounter=1
foreach ii in 32 36 25 35 51 14 11 47 48 {  /* begin loop strata */
scalar jhat=`ii'
di jhat

matrix xe=J(nsimu, 9, 0)
matrix oxed=J(nsimu,1,0)
matrix oxeb=J(nsimu,1,0)

matrix xed=J(nsimu,1,0)
matrix xeb=J(nsimu,1,0)


scalar perpeh=1
scalar jj1=(jhat-1)*8*3+(perpeh-1)*8+1
scalar jj2=jj1+6
matrix a=y[jj1..jj2,1]              /*original data*/
run "$code/analya"
matrix m_mate=mijk[1..8,modeli]           /*estimated parameters*/
scalar modelixe=modeli



scalar perpeh=2
scalar jj1=(jhat-1)*8*3+(perpeh-1)*8+1
scalar jj2=jj1+6
matrix a=y[jj1..jj2,1]              /*original data*/
run "$code/analya"
matrix m_mats=mijk[1..8,modeli]           /*estimated parameters. True values for simulation*/
scalar ms=int(nmv[modeli, 1]+0.5)   				/*nctotal*/
scalar modelixs=modeli


forval inumi = 1/1000{              /*inumi=number of simulations*/

* E
forval i = 1/7{
scalar m=m_mate[`i',1]
matrix a[`i',1]=rpoisson(m)
}
matrix ae=a
scalar m=m_mate[8,1]
matrix xe[`inumi',1]=rpoisson(m)+one7*ae    /*1. Simulated for E */

*Estimate on simulated data
run "$code/analya"
****************
if modeli>0{
matrix xe[`inumi',2]=int(nmv[modeli, 1]+0.5)   				/*nctotal*/
matrix xe[`inumi',3]=nmv[modeli, 2] 				     /*vnctotal*/
*matrix bm000[jhat,2]=smodel[jhat,1]-bm000[jhat,1]  /*scalar vm000=vnctotal-m000*/
} /*modeli*/
****************
* S
forval i = 1/7{
scalar m=m_mats[`i',1]
matrix a[`i',1]=rpoisson(m)
}
matrix as=a
scalar m=m_mats[8,1]
matrix xe[`inumi',4]=rpoisson(m)+one7*as  /*3. Simulated for S */

*Estimate on simulated data
run "$code/analya"
****************
* Direct estimation for the Shining Path
if modeli>0{
matrix xe[`inumi',5]=int(nmv[modeli, 1]+0.5)   				/*nctotal Direct Estimation*/
matrix xe[`inumi',6]=nmv[modeli, 2] 				     /*vnctotal Direct Estimation*/
matrix oxed[`inumi',1]=1
matrix xed[`inumi',1]=xe[`inumi',5]-xe[`inumi',4]



} /*modeli*/
****************
* TRC estimation for the Shining Path
*E+S a
matrix a=ae+as                   /*5. Simulated data of sum of E+S */
*Estimate on simulated data
run "$code/analya"
****************
if modeli>0{
matrix xe[`inumi',7]=int(nmv[modeli, 1]+0.5)- xe[`inumi',2]	 /*nctotal TRC*/
matrix xe[`inumi',8]=nmv[modeli, 2]+ xe[`inumi',3]          /*vnctotal TRC*/

if xe[`inumi',7]-xe[`inumi',4]>0{
matrix oxeb[`inumi',1]=1
matrix xeb[`inumi',1]=xe[`inumi',7]-xe[`inumi',4]
matrix xe[`inumi',9]=xe[`inumi',7]-xe[`inumi',4]
}

} /*modeli*/
****************




} /*simulation*/


*matrix list xe

matrix rmsed=xed'*xed
matrix rmseb=xeb'*xeb
matrix nrmsed=oxed'*oxed
matrix nrmseb=oxeb'*oxeb
scalar rmseds=sqrt(rmsed[1,1]/nrmsed[1,1])
scalar rmsebs=sqrt(rmseb[1,1]/nrmseb[1,1])


matrix vedm=xe[1...,6]
matrix vebm=xe[1...,8]
matrix vedmx=oxed'*vedm/nrmsed[1,1]
matrix vebmx=oxeb'*vebm/nrmseb[1,1]
scalar ved=sqrt(vedmx[1,1])
scalar veb=sqrt(vebmx[1,1])



clear
svmat xe, names(reg)
*summ
summ reg4 if reg4>0,d
scalar msample=r(mean)  
scalar ssample=r(sd)  

summ reg5 if reg5>0,d
scalar mdir=r(mean)  
scalar sdir=r(sd)  

sum reg7 if reg8>0,d
scalar mbasm=r(mean)  
scalar sbasm=r(sd)  

		  
matrix res[icounter,1]=jhat
matrix res[icounter,2]=ms
matrix res[icounter,3]=msample
matrix res[icounter,4]=ssample

matrix res[icounter,5]=mdir
matrix res[icounter,6]=mdir-msample

matrix res[icounter,7]=ved
matrix res[icounter,8]=sdir
matrix res[icounter,9]=rmseds
matrix res[icounter,10]=nrmsed[1,1]
matrix res[icounter,11]=mbasm
matrix res[icounter,12]=mbasm-msample
matrix res[icounter,13]=veb
matrix res[icounter,14]=sbasm
matrix res[icounter,15]=rmsebs
matrix res[icounter,16]=nrmseb[1,1]

scalar icounter=icounter+1			  
} /* end loop strata*/
matrix list res
clear
svmat res, names(v)
ren v1 Stratum		
ren v3 N 
format %7.2f N

ren v5 Estimated 
ren v6 Biasd	
ren v7 Est_SEd	
ren v8 Sample_SEd	
ren v9 RMSEd	
ren v10 Valid_Obsd

ren v11 Estimateb
ren v12 Biasb	
ren v13 Est_SEb
ren v14 Sample_SEb	
ren v15 RMSEb	
ren v16 Valid_Obsb

la var Estimated "Estimate" 
la var Biasd "Bias"	
la var Est_SEd "Estimated S.E."	
la var Sample_SEd "Sample S.D."	
la var RMSEd "Sample RMSE"	
la var Valid_Obsd "Valid Observations"

export excel Stratum N Estimated Biasd Est_SEd Sample_SEd RMSEd Valid_Obsd Estimateb Biasb Est_SEb Sample_SEb RMSEb Valid_Obsb ///
using "$excel/${name}.xls", sheet("res") sheetmodify firstrow(varlabels) missing(".")

capture log close



